The E2 state of FeMoco: Hydride Formation versus Fe Reduction and a Mechanism for H2 Evolution

Abstract The iron‐molybdenum cofactor (FeMoco) is responsible for dinitrogen reduction in Mo nitrogenase. Unlike the resting state, E0, reduced states of FeMoco are much less well characterized. The E2 state has been proposed to contain a hydride but direct spectroscopic evidence is still lacking. The E2 state can, however, relax back the E0 state via a H2 side‐reaction, implying a hydride intermediate prior to H2 formation. This E2→E0 pathway is one of the primary mechanisms for H2 formation under low‐electron flux conditions. In this study we present an exploration of the energy surface of the E2 state. Utilizing both cluster‐continuum and QM/MM calculations, we explore various classes of E2 models: including terminal hydrides, bridging hydrides with a closed or open sulfide‐bridge, as well as models without. Importantly, we find the hemilability of a protonated belt‐sulfide to strongly influence the stability of hydrides. Surprisingly, non‐hydride models are found to be almost equally favorable as hydride models. While the cluster‐continuum calculations suggest multiple possibilities, QM/MM suggests only two models as contenders for the E2 state. These models feature either i) a bridging hydride between Fe2 and Fe6 and an open sulfide‐bridge with terminal SH on Fe6 (E2‐hyd) or ii) a double belt‐sulfide protonated, reduced cofactor without a hydride (E2‐nonhyd). We suggest both models as contenders for the E2 redox state and further calculate a mechanism for H2 evolution. The changes in electronic structure of FeMoco during the proposed redox‐state cycle, E0→E1→E2→E0, are discussed.

. Calculated relative energies in kcal/mol (TPSSh level of theory using the minimal cluster model) and Mulliken spin populations on the metal ions for all 4 brokensymmetry solutions (235, 346, 247, 147-if found) of the E2 models shown in Fig. 3 Figure S1. Calculated relative energies (polarized QM energy from QM/MM calculation) of different E2 models in the ⍺-195 His -Nε(H) protonation state. All models were calculated with 4 different broken-symmetry solutions with the lowest-energy one indicated. See Table S2 for other BS solutions. t3 denotes a proton or a hydride that points to S3A, e.g. S2Bt3 is a protonated S2B where the proton points in the general direction of S3A (instead in direction of S5A) b SH proton points toward homocitrate Figure S2. Calculated relative energies (polarized QM energy from the QM/MM calculation) of different E2 models in the ⍺-195 His -Nδ(H) protonation state. All models were calculated with 4 different broken-symmetry solutions with the lowest-energy one indicated. See Table S3 for additional information.   Figure S3. The electronic structure of the E2-nonhyd model in the BS7-235 solution as interpreted via localized orbital analysis (Pipek-Mezey) and shown as Noodlemanstyle majority/minority spin vectors. Localized orbital isosurfaces (isovalue=0.05) of the minority-spin electrons are shown as insets. Large arrows indicate 5-electron s=5/2 majority-spin vectors, these 5 electrons are well localized. Small arrows indicate minority-spin single s=1/2 electrons that can be more delocalized. Figure S4. The electronic structure of the E2-hyd model in the BS7-247 solution as interpreted via localized orbital analysis (Pipek-Mezey) and shown as Noodleman-style majority/minority spin vectors. Localized orbital isosurfaces (isovalue=0.05) of the minority-spin electrons are shown as insets. Large arrows indicate 5-electron s=5/2 majority-spin vectors, these 5 electrons are well localized. Small arrows indicate minority-spin single s=1/2 electrons that can be more delocalized.  as Noodleman-style majority/minority spin vectors. Localized orbital isosurfaces (0.05 isovalue) of the minority-spin electrons are shown as insets. Large arrows indicate 5electron s=5/2 majority-spin vectors, these 5 electrons are well localized. Small arrows indicate minority-spin single s=1/2 electrons that can be more delocalized. Figure S7. The electronic structure of the E2-nonhyd model in the BS7-247 solution as interpreted via localized orbital analysis (Pipek-Mezey) and shown as Noodlemanstyle majority/minority spin vectors. Localized orbital isosurfaces (0.05 isovalue) of the minority-spin electrons are shown as insets. Large arrows indicate 5-electron s=5/2 majority-spin vectors, these 5 electrons are well localized. Small arrows indicate minority-spin single s=1/2 electrons that can be more delocalized. Figure S8. The electronic structure of the saddlepoint structure (E2-hyd-TS) for H2 formation. The model shown has ⍺-195 His in the Nδ protonation state. Figure S9 shows the effect of increasing the QM-region from 134 atoms to 214 atoms by including the peptide backbone surrounding S3A (featuring weak NH...S hydrogen bonds) on the relative energies of the most stable E2 isomers.

The stability of a protonated carbide model
Two different protonated carbide models (see Figure S10) were tested and compared to the most stable hydride isomer (bH(2,6)-OBS(6) in the BS7-247 state); these were pC(1)-CBS(S2B) in the BS7-346 state and pC(2)-np-CBS in the BS7-346 . Their relative energies reveal that carbide protonation is quite uphill in energy compared to forming a hydride formation. Figure S10. Geometries of protonated carbide models (middle and right) compared to the most stable hydride model (left) as well as polarized QM energies and full QM/MM energies (kcal/mol).

The orientations of hydrogens in E2 isomers
Previous investigations on reduced states of FeMoco, by us and Ryde and coworkers, have found that the direction of protonated sulfides of the FeMoco cluster is affected by a fair amount of the protein environment. The scope of the present study was partially based on previous studies 1,2,3 that investigated proton orientation aspects in various En states. Ryde et al. found that for E1 (using TPSS), S2B protonation is preferred to be oriented towards S3A, the other direction being uphill by ~1.7 kcal/mol, while for S5A the preferred direction is towards S3A, the other direction being uphill by ~5.0 kcal/mol. 3 Protonation of S3A was high in energy, ~13 kcal/mol uphill, with little difference in direction of the proton, towards S5A being ~0.5 kcal/mol uphill.
Our own result for E1 (with TPSSh) agree with this result as we find an energy difference of ~1.4 kcal/mol for different orientations of S2B protonation, however, this is affected by ⍺-195 His protonation scheme.
For E2 models, with the ⍺-195 His -Nε(H) scheme, the S2B protonation is preferred to orient towards S3A (5.7 kcal/mol for E2-bH(2,6)-CBS(S2B)), similar to that found by Ryde and coworkers, but with the ⍺-195 His -Nδ(H) scheme the S5A direction becomes favoured instead, with the S3A direction 1.3 kcal/mol less favorable for E2-bH(2,6)-CBS(S2B). This has to do with favourable hydrogen bonds to ⍺-195 His , as shown in Figure S11. Similarly, the S5A protonation is also found to be preferred towards S3A, ~6.3 kcal/mol uphill for E4 models. This is due to the hydrogen bonds that are formed with ⍺-359 Arg and ⍺-96 Arg , though the hydrogen bond with ⍺-96 Arg is more important in this context, see Figure S12. Protonating S3A is always ~10-20 kcal/mol uphill in energy due to the steric effects of the nearby protein backbone as shown in Figure S13.
As hydrides are also present in the E2 models, we also tried different directions for hydrides, terminal hydrides that are trans to the interstitial carbide vs. endo hydrides (with an acute angle with respect to the central carbide while not bridging two Fe ions), bridging hydrides that point towards different sides of the dihedral plane formed by the central carbide, two binding Fe and µ 2 -S. Terminal hydrides trans to the carbide were generally favoured, from 1.2-8.2 kcal/mol depending on what Fe forms a hydride and what sulfide gets protonated, though terminal hydrides are generally unfavored. With the bridging hydrides, as in the case of protonated belt-sulfides, the direction preferred was dependent on the ⍺-195 His protonation scheme. For the ⍺-195 His -Nε(H) scheme the different directions are approximately equal in energy. For the ⍺-195 His -Nδ(H) protonation scheme the hydride is preferred to be oriented towards S5A by ~4.3 kcal/mol (E2-bH(2,6)-CBS(S2B)). Lastly, we investigated rotations of the open protonated belt sulfides and found a preference for orienting towards the ⍺-195 His amino acid residue for ⍺-195 His -Nε(H) scheme and towards the homocitrate ligand for the ⍺-195 His -Nε(H) scheme in the case of the lowest hydride structure (E2-hyd). As the dihedral angle of the hydrides for the OBS models is close to planar, there is no preference involved for the hydrides.
While we find these hydrogen orientations to be non-negligible and we tried to probe the most relevant orientation in each case, it was not possible to try every conceivable permutation of protons/hydride on the cofactor. In our experience this orientation dependence has a small effect on the electronic and molecular structure of the cofactor but rather reflects the H-bonding network surrounding the cofactor. Future work may involve a dynamics-based approach to further study this hydrogen orientation dependence.

The effect of basis set and density functional choice on stability of models
In order to test the sensitivity of our computational protocol with respect to basis set size and functional choice we recalculated the energies (full geometry optimizations with the cluster-continuum model) with either a non-hybrid functional TPSS or the hybrid-functional TPSSh (results were calculated without ZORA ) and either a double-zeta basis set def2-SVP or a triple-zeta basis set def2-TZVP (should be reasonably close to the basis set limit). The results are seen in Table S8. They reveal that there is a basis set dependency present that can affect the relative energies by a few kcal/mol that in turns affects slightly the energetic ordering of isomers. A larger effect is seen by changing the functional from TPSSh to TPSS, which reveal that together with basis set effects, different computational protocols can lead to complete disagreement about what the most stable isomer for E2 is. As found in our previous work (Thorhallsson et al. Chem. Sci. 2019, 10, 11110-11124), the TPSSh functional (using a triple-zeta basis set) leads to much better agreement with the high-resolution crystal structure than TPSS. Table S8. Calculated relative energies (kcal/mol) various hydride and non-hydride E2 models using the cluster-continuum models with different functionals (ZORA was not included in this comparison).

Localized orbital populations
Tables S9-S15 show populations of the Pipek-Mezey localized orbitals for the indicated BS-state models.      The lowest energy models of each structural class in Figure 3 (bottom row) were chosen (with the exeption of bH(2,6)-CBS(S2B)) and ~28 BS solutions were calculated (full geometry optimization with ZORA-TPSSh in a CPCM continuum). BS solutions 123,124,134,234 and 567 were skipped as they were found to always involve an unfavorable Mo(III) Hund configuration resulting in an unfavorable state energy of >20 kcal/mol (see bH(4,5)-OBS(5)-134 in Table S16 for an example). Energies and Mulliken spin populations of all states are shown in Tables S16-S19.  Table S17. Calculated relative energies in kcal/mol (TPSSh level of theory using the minimal cluster model) and Mulliken spin populations on the metal ions for all brokensymmetry solutions of the E2 bH(4,5)-OBS(5) models shown in Fig. 3 Table S18. Calculated relative energies in kcal/mol (TPSSh level of theory using the minimal cluster model) and Mulliken spin populations on the metal ions for all brokensymmetry solutions of the E2 noH-CBS(S2B,S5A) models shown in Fig. 3 Table S19. Calculated relative energies in kcal/mol (TPSSh level of theory using the minimal cluster model) and Mulliken spin populations on the metal ions for all brokensymmetry solutions of the E2 tH(5)-CBS(S2B) models shown in Fig. 3